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SUMMARY 

The  parabolic  equation  is  a  full-wave  approximation  to  the  Helmholtz  wave  equation  that, 
unlike  existing  ray  tracing  and  coupled  mode  techniques,  readily  accommodates 
(wo-dimensionally  inhomogeneous  atmospheres.  This  report  discusses  the  implementation 
of  a  split-step  Fourier  Transform  solution  to  the  parabolic  equation  and  includes  examples 
of  coverage  diagrams  calculated  by  this  method. 
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GLOSSARY  OF  METEOROLOGICAL  TERMS 

Horizontal  movement  of  an  air  mass  mainly  caused  by  wind 
movement  around  high  or  low  pressure  systems.  Super-refraction 
due  to  advection  occurs  primarily  on  over-sea  paths  as  an  air  mass 
that  has  become  warm  and  dry  over  land  is  carried  out  to  sea 
where  it  overlies  a  layer  of  cool  moist  air,  establishing  both  a 
temperature  inversion  and  an  increased  humidity  lapse  rate  in  the 
boundary  region.  The  effect  can  operate  in  reverse  as  when  cool 
moist  air  flows  from  the  sea  onto  warm  dry  land. 

Boundary  region  formed  between  warm  dry  air  stream  above  a 
cool  stream  caused  when  air  heated  over  the  land  surface  rises  and 
is  replaced  by  cooler,  denser  air  from  the  sea  so  that  both  a  sea 
breeze  near  the  surface  and  a  strong  off-shore  movement  of  dry  air 
is  created.  Generally  extending  above  100  km,  coastal  fronts  are 
common  along  the  western,  southern  and  eastern  Australian 
coasts. 

Effect  common  in  dear-sky  conditions  over  land  paths  resulting 
from  a  temperature  inversion  created  when  air  near  the  ground  is 
cooled  considerably  whilst  the  warmer  upper  air  (heated  via  solar 
radiation  during  the  day)  remains  warm.  Cloud  cover  reflects 
radiation  back  to  the  ground,  thereby  avoiding  such  cooling. 
Nocturnal  cooling  does  not  occur  over  sea  paths  as  diurnal  ocean 
temperatures  do  not  vary  by  much. 

Effect  due  to  an  air  mass  which  descends  and  hence  is  compressed 
and  heated  adiabatkally.  The  descending  air  is  very  dry  since  it 
originated  at  a  high  altitude  and  so  an  increased  humidity  lapse 
and  possibly  a  temperature  inversion  occurs.  Subsidence 
inversions  are  more  frequent  and  occur  at  tower  altitudes  in  winter 
than  in  summer. 
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CCIR 

ERL 

International  Radio  Consultative  Committee 

Electronics  Research  Laboratory 

FFT 

Fast  Fourier  Transform 

IREPS 

Integrated  Refractive  Effects  Prediction  System 

ODE 

Ordinary  Differential  Equation 

VHF 

Very  High  Frequency 
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1  INTRODUCTION 

r' 

Under  the  influence  of  synoptic  processes  such  as  advent**,  subsidence,  coastal  fronts,  or  nocturnal 
cooling,  stratification  of  the  troposphere  is  the  focus  of  re fr activity  layering  can  occur,  if  the  resulting 
refraction  of  propagating  electromagnetic  waves  is  sufficiently  peat  to  cause  a  ray  curvature  that 
exceeds  the  curvature  of  the  earth's  surface,  then  the  propagating  wave  can  be  channelled  in  a  duct. 
At  VHF  and  higher  frequencies,  surface  and  elevated  tropospheric  ducts  can  cause  extended 
propagation  well  beyond  the  radio  horizon  with  little  attenuation  relative  to  free  space.  The  received 
radio  djm.k  can  then  be  markedly  stronger  (or  weaker)  than  would  be  opected  for  s  standard 
atmosphere.  The  anomalous  propagation  of  energy  via  a  dud  can  thereby  lead  to  the  volatile 
problem  of  inter-system  interference,  whilst  radar  detection  ranges  can  be  enhanced  for  targets 
located  within  the  duct.  However,  ducts  can  also  decrease  radio  (radar)  operational  ranges  by 
causing  "holes"  in  signal  coverage  if  a  receiver  (target)  is  located  outside  the  duct.  Reliable  prediction 
of  such  anomalous  propagation  effects  can  therefore  enhance  the  effective  deployment  of 
communications,  surveillance  and  electronic  warfare  systems. 

Propagation  prediction  models  tend  to  fall  into  two  categories:  those  methods  based  on  optical  ray 
tracing  techniques;  and  those  relying  on  mode  theory.  Although  qualitative  predictions  can  be  made 
using  ray  theory,  problems  associated  with  the  focusing  of  rays  and  identification  of  ray  families 
render  these  methods  unsuitable  for  modelling  diffraction  effects  in  non-standard  atmospheres. 
Mode  theory,  which  involves  a  full-wave  solution  10  Maxwell’s  equations,  can  be  applied  to  ducting 
problems  but  becomes  intractable  when  complex  refractive  structures  arc  involved.  Computing  cost 
lends  to  limit  convergent  solutions  for  modal  analyse,  to  ample  refractive  structures.  Neither  of  these 
methods,  nor  existing  hybrid  combinations  of  these,  can  be  readily  applied  lo  the  two- dimensionally 
inhomogeneous  atmospheres  of  range-dependent  problems. 

Since  the  wave  equations  governing  realistic  physical  representation  of  the  troposphere  permit  a 
closed  form  solution  only  in  very  simple  cases,  appi  oxiroatioos  to  these  equations  have  been 
introduced  that,  with  the  availability  of  modem  uumenca-  'echniqucs  and  fast  computers,  facilitate 
solution  of  complicated  propagation  problems.  This  re;  on  demonstrates  the  application  of  an 
alternative  full-wave  model  that  is  based  upon  such  an  approximation  to  the  wave  equation  and  can 
be  applied  to  complicated  range-dependent  tropospheric  propagation  problems.  The  parabolic 
equation  was  first  obtained  for  propagation  over  a  spherical  earth  in  a  vertically  inhomogeneous, 
horizontally  homogeneous  atmosphere  by  (Ueontovich  and  Fock,  1946)  and  has  been  extended  by 
(Ko  rt  of,  1964)  for  the  case  of  a  two-dimension  ally  inhomogeneous  atmosphere  lo  give  a  parabolic 
partial  differential  equation  for  the  field  amplitude.  The  model  presented  in  this  report  uses  a  split 
step  Fourier  transform  technique  which  takes  advantage  of  the  computational  speed  of  the  Fast 
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Found  Trattfom  ttg  ,  am  and  is  numerically  stable.  Thi*  solution  method  was  developed  in  the 
underwater  acoustics  arena  by  (Hardin  and  Tappen.  1973)  to  solve  a  paraboLc  equation  of  identical 
form  to  'he  parabolic  equation  for  tropospheric  propagation.  The  use  of  this  technique  affords  a  full- 
wave  solution  without  the  large  computational  overheads  of  previous  fuff-wave  methods,  ft  is  this 
property  of  the  method  which  makes  k  most  attractive. 

This  r*pott  begins  with  a  brief  description  of  the  wave  equation  governing  two-dimensional 
tropospheric  propagation  and  indicates  how  propagation  problems  can  be  represented  by  a  partial 
differential  equation  vhh  camples  coefficients  After  a  tfr^utrino  of  the  applicability  of  existing 
prediction  techniques  for  ducting  environments,  it  is  shows  how  the  tropospheric  wave  equation 
reduces  to  a  parabolic  equation  that  facilitates  ready  solution  This  approximation  to  the  wave 
equation  has  the  significant  advantage  that  it  invites  'marching’  type  solutions  since  the  partial 
different  ui  equation  is  now  orty  fir  f  order  in  range 


The  third  section,  "Split  step  Fourier  Solution",  is  devoted  to  a  description  of  the  solution  method 
chose-  to  evaluate  the  parabolic  equation  of  section  2.  The  split  step  solution  is  a  'marching"  type 
solution  since  the  solution  can  be  evaluated  at  range  x+$x  (4x  is  sufficiently  small)  by  using  the 
evaluated  solution  at  range  mu  initial  value  far  calculations  Providing  that  appropriate  boundary 
conditions  are  satisfied  the  solution  can  in  this  way  be  stepped  out  to  the  required  range-  The  entire 
range-  and  height -dependent  field  is  calculated  as  the  solution  is  marched  forward  in  this  manner. 
The  solution  involves  the  discretisation  of  the  electromagnetic  field  with  respect  to  height  at  each 
range  step  and  the  use  of  the  Fast  Fourier  Transform.  In  this  section,  solution  stability  is  considered 
and  an  indication  of  the  errors  inherent  in  the  technique  is  also  given. 

Since  a  marching  solution  can  be  applied  to  the  parabolic  equation,  a  maximum  range  boundary 
condition  need  not  be  specified,  thus  simplifying  computation.  Tropospheric  propagation  problems 
can  therefore  be  specified  as  initial  boundary  value  problems  on  an  open  domain.  Section  4  deals 
with  the  specification  of  this  initial  starting  field  and  of  the  surface  and  maximum  altitude  boundary 
conditions.  There  arc  however  some  signal  processing  aspects  of  implementation  that  need  to  be 
addressed,  and  included  in  this  section  is  a  discussion  of  these.  The  following  section  contains  an 
analysis  of  solution  sensitivity  to  perturbations  in  initial  field  and  to  small-scale  discontinuities  in 
refractive  index  profiles. 

Within  a  simplified  ray  tracing  context,  section  6  illustrates  the  effects  that  variations  in  refractive 
structures  have  upon  tropospheric  propagation.  The  main  purpose  of  this  section  is  to  introduce  the 
categories  of  ducting  profile*  that  art  modelled  in  section  7. 
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The  report  concludes  with  representative  examples  of  coverage  diagrams  generated  using  the  split 
step  solution  to  the  parabolic  wave  equation.  These  diagrams  are  im.ructive  as  examples  which 
illustrate  the  complicated  nature  of  coverage  patterns  resulting  from  common  refractive  index 
structures.  They  were  generated  on  a  COMPAQ  386/20e  personal  computer  with  an  Eighteen  Eight 
Labs  fast  parallel  processing  board.  The  program  used  to  generate  the  diagrams  in  this  report  is 
currently  under  development  in  Communicalions  Division,  Electronics  Research  Laboratory  (ERL) 
at  Defence  Science  and  Technology  Organisation.  It  accepts  environmental  data  for  the  propagation 
paths  of  interest  in  the  form  of  refractivity  profiles  (specified  either  analytically  or  by  measurement 
data)  and  generates  two-dimensional  Height-Range  coloured  intensity  plots  for  the  area  of  maximum 
height  and  range  specified  by  the  operator.  The  program  has  also  been  used  to  generate  two-way 
path  toss  diagrams  for  radar  applications  (Slingsby,  1990). 

In  addition  to  the  large  computational  saving  over  other  full-wave  methods,  a  major  advantage  of  the 
parabolic  equation  technique  outlined  in  this  report  is  that  it  affords  progressive  calculator  of  the 
propagated  field  from  one  range  location  to  the  next.  Hence  the  ERL  developmental  software  has 
potential  application  as  an  interactive  research  tool  by  which  an  operator  would  be  able  to  create 
desired  propagation  patterns  through  interactive  manipulation  of  the  input  environmental 
parameters.  This  is  an  advantage  which  other  applicable  techniques  that  have  been  developed  to  dale 
do  not  have.  For  those,  the  environmental  input  information  needs  to  be  defined  before  computation 
commences.  Thus  given  a  means  of  remotely  sensing  the  atmospheric  refractivity  data,  the  ERL 
developmental  software  could  be  adapted  to  provide  a  real-time  updating  display  illustrating  the 
effects  of  the  transient  weather  conditions  on  signal  propagation. 


2  TROPOSPHERIC  WAVE  PROPAGATION  PROBLEMS 

This  section  formulates  the  tropospheric  wave  propagation  problem  and,  after  a  brief  discussion  of 
prediction  techniques  currently  available  for  modelling  ducting  effects,  introduces  tlie  parabolic 
approximation  to  the  elliptic  wave  equation.  Consideration  of  the  significance  of  this  ap  proximal  ion 
for  tropospheric  propagation  problems  indicates  that  for  realistic  ducting  envir  mments  the 
assumptions  inherent  in  this  parabolic  approximation  are  valid.  Whereas  it  usual  nurocri  a>  schemes 
the  elliptic  wave  equation  requires  soiutioa  of  a  large  system  of  simultaneous  equation.'  i  a  large 
number  of  unknowns  with  boundary  conditions  being  specified  on  a  dosed  domain,  t'.c  parabolic 
equation  can  instead  be  solved  by  a  ’marching’  technique  on  an  open  domain;  so  numerical  solution 
is  easier  to  compute. 
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2.1  Governing  wave  equation 

The  Helmholtz  wave  equation  for  the  electric  or  magnetic  field,  ft,  of  an  isotropic  time- 
harmonic  point  source  is  given  by 

V1  *  +  k1  n1  ♦  =  0  (1) 

where  V  is  the  Lapladan  operator,  k=2w  /X  is  the  free  space  -wavenumber  and 
nmJ(e  -ja)/w€0is  the  index  of  refraction,  with  £  0  denoting  free  space  permittivity  and  E ,  a 
representing  the  absolute  permittivity  and  conductivity  respectively  of  the  propagating 
medium. 

In  order  to  examine  the  spherical  spreading  behaviour  from  e  point  source  of  electromagnetic 
waves  along  a  spherical  earth  it  is  convenient  to  adopt  the  spherical  coordinate  system  (r,0  fp ) 
with  the  origin  at  the  centre  of  the  earth  and  with  the  polar  axis  drawn  through  the  source 
(0  >0  at  the  source).  The  first  assumption  to  be  made  in  this  study  is  that  of  azimuthal 
symmetry.  This  is  introduced  for  two  reasons:  firstly,  it  is  rare  that  sufficient  environmental 
information  is  available  to  warrant  the  computational  overhead  for  a  full  three-dimensional 
solution  to  the  wave  equation;  and  secondly,  it  is  assumed  that  the  azimuthal  variations  of 
refractive  index  will  be  sufficiently  slow  to  justify  a  two-dimensicnal  approach.  The  field,  ft , 
and  refractive  index,  n,  will  in  this  case  be  functions  of  (rj0 )  t  nly.  The  tropospheric  wave 
propagation  problem  can  therefore  be  viewed  in  terms  of  a  closed-domain  boundary  value 
problem,  as  indicated  by  figure  1. 


Figure  1  Tropospheric  boundary  value  problem 
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The  assumption  of  azimuthal  symmetry  means  that  (  d/d$  )s  0  in  Maxwell's  equations.  For 
an  electric  dipole,  then,  only  field  components  Ef,  Eq  and  will  be  non-zero;  and  for  a 
magnetic  dipole,  only  Hf,  Hq  and  will  exist.  In  the  case  of  an  electric  dipole,  Maxwell's 
equations  reduce  to 

7fa-,^>-ir)+^o^=0 
r-iella(H«^>>=^Er 
•7  ;,<&*>  =  **** 

and  combining  these  three  equations  to  eliminate  field  components  Ef  and  Eq 

-  <  |  -  f-(rH , )  1  -  ;  £  [— ~  pI(ha  sine)I  }  "  “2RnHA  =  0 

r  dr  c  dr  <p  r  do  c  sin©  cO  9  0  <p 

<?  d 

Expanding  (2),  using  the  subscripts  r  and  0  to  denote  —  and  ^  ^ 
dropping  the  <p  subscript  from  the  scalar  H^: 


(2) 


g-g  respectively,  and 

2tlr  H  _  Hcote  H# cote  _  H.e,  He,_  H  cote  c,  _  H,cfl  f  _  Q 

,r  r  ”  r1  r*  r‘  c  rc  r*c  r*  r‘c  v  ' 


since 

and 


c  =  c  n 
o 

2  i  2 

u  u  c  =  k 
*0  o 


Equation  (3)  can  be  represented  in  the  form 


VZH  ♦  k2n*H .  +  ~  X  (  V  X  H,)  =  0  (■») 

p  <p  C  p 

22  Existing  prediction  techniques 

There  exists  a  large  body  of  literature  devoted  to  the  theory  of  ray  tracing  and  coupled  mode 
techniques,  and  so  these  will  receive  only  a  cursory  treatment  here.  Hybrid  combinations  of 
ray  tracing  and  mode  methods  also  exist;  perhaps  the  most  widely  used  of  these  is  the 
Integrated  Refractive  Effects  Prediction  System  (IREPS)  developed  by  the  Naval  Ocean 
Systems  Centre,  San  Diego.  However,  [REPS  does  not  overcome  many  of  the  limitations 
inherent  in  ray  tracing  and  mode  methods  used  for  ducting  applications. 

The  main  purpose  of  this  section  is  to  highlight  the  deficiencies  of  such  techniques  thereby 
emphasising  the  advantage  to  be  gained  through  implementation  of  the  alternative  parabolic 
equation  method. 
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2.2.1  Geometrical  optics 

The  general  geometrical  theory  can  be  found  in  a  number  of  works.  See,  for  example, 
(Budden,  1961), 

Assuming  that  the  refractive  index  of  the  medium  is  slowly  varying  in  space,  then  an 
approximate  solution  to  the  wave  equation  (1)  is  sought  in  the  form 

«(r)  =  A(r)e^r>  (5) 

where  A(r)  and  S( r)  are  slowly  varying  functions  of  position  vector  r;  the  function  S, 
known  as  the  icoual  function,  defines  surfaces  of  constant  phase.  For  this  solution  to 
satisfy  the  reduced  wave  equation  it  is  necessary  that 

!ys|*  =  n*(r)  (6) 

and 

A(r)  v  S(r)  +  2  \7S(r)  A(r)  =  0  (7) 

hold.  The  surfaces  S(r)  =  constant  are  the  wave  fronts  and  so  the  ray  paths  along  which 
the  electromagnetic  waves  propagate  are  given  by  the  trajectories  orthogonal  to  these 
surfaces  (ic  |  \7s| 1  ).  Integration  of  equation  6  (known  as  the  iconal  equation)  and  of 
equation  7  give  ray  trajectories  and  wave  amplitude. 

If  the  atmospheric  refractive  index  can  be  assumed  to  be  spherically  stratified  with 
respect  to  the  surface  of  the  earth  (it  is  horizontally  homogeneous,  varying  with  altitude 
only),  then  equation  (6)  is  equivalent  to  Snell’s  law  and 


n(h)  sin  9  ■  S 


(8) 


where  8  is  the  angle  that  a  ray  makes  with  the  vertical,  S  is  a  constant  for  each  ray,  and 
h  refers  to  altitude  measured  from  the  earth’s  surface.  The  radius  of  curvature  of  a  ray, 
p,  is  then  given  by  (Livingston,  1970) 


P 


•n 

&  sin  6 


(9) 


There  are  several  limitations  to  ray  tracing  techniques.  The  first  of  these  relates  to  an 
assumption  that  the  fractional  change  in  spacing  between  neighbouring  rays  (initially 
parallel)  is  small  with  respect  to  a  wavelength.  In  a  ducting  environment,  where  rays 
are  ’trapped’  and  are  strongly  refracted  or  reflected,  this  requirement  will  be  violated. 


f 


6 


UNCLASSIFIED 


UNCLASSIFIED 


ERL-0531-RR 


| 

Another  limitation  of  ray  tracing  techniques  involves  the  identification  of  ray  families. 
A  ray  is  assumed  to  be  a  typical  member  of  a  ’family'  that  consists  of  a  set  of  rays  with 
similar  initial  conditions  and  which  have  followed  similar  paths.  In  order  to  calculate 
the  field  at  a  given  point,  the  contribution  of  each  ray  passing  through  that  point  must 
be  determined,  and  therefore  the  ray  families  must  be  identified.  Errors  occur  when 
rays  are  incorrectly  identified,  is  which  case  their  contribution  to  the  total  field  will  be 
inaccurately  calculated. 

At  caustics  (surfaces  occurring  where  adjacent  rays  from  the  same  family  intersect),  the 
above  ray  equations  wilt  have  multiple  solutions;  and  at  extended  ranges,  as  more 
caustics  are  generated,  ray  theory  tends  to  break  down  completely.  Asymptotic 
approximations  for  the  vicinity  of  caustics  tend  *o  be  complex  and  to  date  no  simple 
method  has  been  developed  that  would  encourage  ready  implementation  of  such 
approximations. 

Frequency  does  not  enter  into  ray  tracing  treatments  of  propagation  in  ducting 
environments  and  so  errors  are  introduced  due  to  the  dependence  on  frequency  in 
determining  whether  a  ray  will  be  trapped  by  a  given  duct  structure.  In  addition,  it  is  a 
requirement  of  fay  theory  that  the  path-length  difference  between  direct  and  reflected 
rays  be  at  least  one-quarter  wavelength  (Fishback,  1951).  Ray  theory  methods  are 
therefore  limited  to  regions  within  the  radio  horizon. 

Ray  tracing  methods  follow  very  simple  geometry  and  yield  a  convenient  qualitative 
picture  of  ducting  phenomena.  However  the  tendency  of  ray  theory  to  reflect  all  or 
nothing  of  the  rays  from  ducting  layers  is  another  limits  ion  of  these  techniques,  as 
appreciable  ducting  is  only  indicated  when  both  transmitting  and  receiving  antennas  are 
within  the  duct.  By  ray  theory,  a  wave  originating  exterior  to  a  duct  generally  cannot  be 
trapped  within  the  duct,  and  there  is  no  allowance  for  leakage  from  the  duct  by  trapped 
waves.  Partial  reflection  from  ducting  layers  has  not  yet  been  adequately  addressed. 
Similarly  lateral  changes  in  refractive  index  have  been  considered  by  a  small  number  of 
authors  although  results  have  been  of  a  qualitative,  rather  than  quantitative  nature. 

2JJ  Coupled  Mode  Techniques 

Since  ray  tracing  theory  does  not  accommodate  non-standard  atmospheres  well,  in 
dealing  with  dueling  environments  a  full-wave  solution  to  Maxwell's  equations  is  often 
sought.  The  waveguide  model  for  electromagnetic  propagation  in  a  tropospheric  duct 
treats  such  propagation  as  being  analogous  lo  propagation  in  a  leaky  waveguide.  For 
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comprehensive  treatment  of  such  methods  refer  to  (Booker  and  Walkinshaw,  1946), 
(Kerr,  1951)  and  (Budden,  1961). 

In  the  tropospheric  wave  problem,  as  formulated  in  section  2.1,  the  surface  of  the  earth 
acts  as  a  curved  boundary  from  which  electromagnetic  waves  may  be  reflected  and 
around  which  these  waves  may  be  diffracted.  If  the  refractive  index,  n(h),  is  assumed  to 
be  spherically  stratified  with  respect  to  the  surface  of  a  smooth  earth,  then  in  terms  of 
an  earth  flattened  geometry,  earth  curvature  is  included  in  a  modified  index  of 
refraction  as  defined  by 

m(h)  =  n(h)  +  |  (10) 

where  h  denotes  altitude  measured  from  the  earth’s  surface,  and  a  is  the  radius  of  the 
earth.  The  field  of  a  source  in  a  stratified  atmosphere  can  then  be  represented  by  a 
sum  of  leaky  waveguide  modes.  For  a  dipole  source,  the  field  will  be  proportional  to 


J-  £  e4kSRFt(h)F  (h)e  (11) 

“  (Budden,  1961) 

where  R  is  the  horizontal  range,  k  is  the  free  space  wavenumber,  e  represents  the 
excitation  of  the  mode,  and  F,(h),  Ff(h)  determine  the  attenuation  of  the  mode  with 
height  (normalised  to  unity  at  the  earth’s  surface)  for  transmiller  and  receiver 
respectively.  The  terms  F(h)  can  be  expressed  in  terms  of  the  height-gain  function 
g(h,S):  F(h)  -  g(h5)/g(0,S),  where  g(h,S)  is  the  outgoing  wave  solution  of 

+  V  [m>  (h)-  S'  ]g(h,S)  -  0  (12) 

Since  S  is  complex,  the  term  exp  ltSR  in  equation  (11)  determines  the  attenuation  of  the 
mode  with  range. 

In  order  to  determine  the  allowed  modes  it  is  necessary  to  apply  boundary  conditions. 
Firstly,  the  radiation  boundary  condition  requires  that  waves  are  only  outgoing.  If  the 
changes  in  refractivity  gradient  associated  with  ducting  layers  are  represented  in  a 
vertical  modified  refractive  index  profile  by  discontinuities  in  a  linear  segmented 
profile,  then  this  radiation  condition  puts  a  constraint  on  the  height-gain  function  at 
heights  above  the  discontinuity  of  a  ducting  layer.  Secondly,  the  boundary  condition  at 
the  earth’s  surface  must  be  satisfied.  This  constrains  the  height-gain  function  in  the 
region  below  a  ducting  layer.  Finally,  the  height-gain  functions  and  their  derivatives 
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above  and  below  the  discontinuities  at  the  boundaries  of  a  ducting  layer,  are  required  to 
be  continuous  across  those  discontinuities.  This  last  boundary  condition  leads  to  the 
mode  equation,  the  details  of  which  can  be  found  in  most  papers  dealing  with  coupled 
mode  theory.  The  reader  is  referred  to  (Craig,  1985)  for  a  useful  interpretation  of  the 
roots  of  the  mode  equation  in  terms  of  the  familiar  ray  theory  picture. 

Locating  the  roots  of  the  modal  equation  is  a  complex  task  and  much  research  has  been 
channelled  towards  this.  Various  numerical  solution  techniques  have  been  investigated 
in  respect  to  this  problem  but  these  will  not  be  discussed  here.  The  reader  is  referred 
to  (Rotheram,  1983)  for  a  survey  of  existing  solution  methods. 

Perhaps  the  most  significant  disadvantage  associated  with  modal  techniques  is  the 
complexity  of  solutions.  The  number  of  modes  which  need  to  be  considered  for  ducting 
problems  is  proportional  to  frequency  and  to  the  thickness  of  the  duct.  Modal  methods 
can  therefore  become  intractable  at  high  frequencies  and  when  complicated  duct 
profiles  are  involved. 

Experimental  verification  of  waveguide  models  has  been  largely  qualitative. 
Undoubtedly  this  reflects  the  difficulty  of  adequately  monitoring  refractivity  structures 
during  experimental  studies,  but  it  also  indicates  some  inadequacy  of  the  modal  ducting 
model  Observational  data  indicate  received  fields  appreciably  stronger  than  can  be 
accounted  for  by  the  simple  homogeneous  duct  structures  considered  in  most  mode 
analyses.  A  number  of  authors  (see,  for  example.  Barton,  1973)  have  suggested  the 
presence  of  ducting  layers  that  vary  with  range  as  a  cause  of  such  observed  phenomena. 
Although  there  have  been  some  attempts  to  address  problems  involving  refractive 
structures  which  vary  with  range  as  well  as  with  height  (refer  to  Wait,  1980,  for 
example),  the  results  have  again  been  qualitative  in  nature.  Computation  is  also 
considerably  increased  in  such  approaches  over  that  Tor  the  homogeneous  ducting 
problems.  Thus  most  implementations  restrict  the  refractive  index  to  be  a  simple 
analytical  function  of  height  alone,  and  therefore  have  application  only  to  range- 
independent  problems.  The  most  common  problem  of  this  type  involves  low-level 
evaporation  ducts  over  the  ocean.  For  problems  involving  surface  and  elevated  ducts, 
however,  modal  models  tend  not  to  be  applicable. 

2.2.3  I  REPS 

The  Integrated  Refractive  Effects  Prediction  System  (Hitncy  and  Richter,  1976)  has 
been  used  since  1978  by  the  United  States  Navy  and  other  organisations  to  give  an 
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operational  propagation  assessment  capability  for  naval  surveillance,  communications, 
electronic  warfare  and  weapon  guidance  systems.  It  is  designed  to  deal  with 
evaporative  and  surface-based  ducting,  and  so  does  not  model  well  the  propagation 
effects  of  elevated  ducts. 

The  products  IREPS  produces  are  based  upon  a  combination  of  ray-optics  and  a 
simplified  full-wave  single-mode  solution,  with  semi-empirical  formulations  based  on 
measured  data  and  interpolations  to  smooth  transition*  between  the  various  models. 
Polarisation  effects  are  neglected  by  IREPS  on  over-the-horizon  paths  which,  like  the 
methods  discussed  in  sections  2.2.1  and  2.2-2,  assumes  range-independence  of  the 
ducting  structures.  An  example  of  a  coverage  diagram  generated  by  IREPS  is  included 
at  figure  2,  where  the  signal  coverage  for  a  homogeneous  duct  between  SO  m  and  600  m 
has  been  modelled.  (The  duct  modelled  here  is  of  identical  strength  and  extent  to  that 
used  for  figure  8,  page  31,  in  which  the  parabolic  equation  method  was  employed  to 
generate  signal  coverage  from  a  gaussian  beam  antenna.) 

For  airborne  applications  particularly,  elevated  ducts  assume  a  greater  signilicance  in 
the  prediction  of  propagated  signals.  For  these  ducts  it  is  no  longer  adequate  to  assume 
horizontal  homogeneity  of  the  inversion  layers.  Real-time  modal  calculations  are 
prohibitive  due  to  the  computational  intensity  involved  in  such  analysis;  whilst  the 
approach  adopted  by  IREPS  of 'template  matching"  -  whereby  precalculated  results  arc 
scaled  to  the  frequency  and  duct  height  of  interest  -  would  be  very  difficult  to  extend  for 
the  case  of  elevated  ducting  layers  (since  the  single-mode  approximation  currently  used 
for  beyond-tbe-horizon  ranges  would  no  longer  be  valid  for  these  elevated  ducts,  for 
which  multi-mode  propagation  is  important).  For  such  reasons,  the  approach  of  the 
parabolic  equation  and  split-step  Fourier  solution  becomes  very  attractive. 
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Figure  2  IREPS  coverage  diagram  far  600  m  homogeneous  duct 
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23  Parabolic  Equation 


By  making  the  substitution 


u  (tB) _ u(r,e)cik*e 

*V(r#)  r  C(rj9)./sin0 


(13) 


as  suggested  by  (Ko  et  al,  1984)t  equation  (3)  becomes 


,  3E-  C„x  2k  4 

u  -  U  {  — T  *  -  >  ♦  U  ( — “ 
r  C  C  Or 


cot©  Ui|  .  3  »  co 

!*•  ra  C  2r*  4r 


I  k  a  cot0 


3  col©  C, 


3Ct 
r  x  C 

x  .  ili*-  See  *  2t-' )2  -  *  kVl  =  0  (U) 

»■  r  C  r*C  »■  r  •  ** 


♦  f 


The  effect  of  the  substitution  of  equation  (12)  is  to  remove  the  rapid  phase  variation  of  the 
field  H^.  The  id  (o'  -  1)  term  in  (13)  arises  from  the  approximation  Id  a*  /  H  K  Id  since,  in 
the  troposphere,  altitude  above  the  earth  will  be  small  compared  to  the  radius  of  the  earth,  a. 
Equation  (14)  can  then  be  reduced  to 

un  +  2ikux  +  ld(n1.l+-%)u  =  0  (15) 

when  the  following  approximations  are  assumed:- 


ka'te’ 

*  nee 

(3) 

e> & 

ka 

(c) 

*e< 

f  - 

ka 

100 

(*») 

€  ,2k3  ,3 

€  '  a  ' 

(d) 

and  the  following  earth-flattening  substitutions  are  incorporated:  - 


z  =  r  -  a  (16) 

x  «  a8  (17) 

o(x,z)  =  n(r£)  +  ^  (18) 

Equation  (IS)  is  known  as  the  parabolic  equation  since  it  is  second-order  in  altitude,  z,  and 

first-order  is  range,  x.  The  advantage  that  an  equation  of  (his  form  has  over  that  of  equation 
(3)  is  that  solution  to  the  parabolic  equation  does  not  require  that  the  maximum  range 
boundary  be  specified.  Subsequently,  the  parabolic  equation  invites  a  ‘marching'  solution  in 
which  the  equation  is  solved  al  initial  range  x*0  and  is  (ben  solved  for  small  steps  out  is  range, 
using  the  solution  at  the  previous  step  as  the  starting  field  for  the  present  calculation.  In  this 
way  the  solution  can  be  calculated  out  to  any  range,  providing  that  initial  field,  upper  and 
surface  boundary  conditions  are  satisfied.  Such  a  solution  is  inherently  easier  to  compute  than 
one  which  requires  a  large  number  of  unknowns  with  boundary  conditions  being  specified  on  a 
closed  domain,  as  is  the  case  for  the  full  elliptic  wave  equation.  The  tropospheric  wave 
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propagation  problem  can  now  be  considered  in  terms  of  the  open-domain  boundary  value 
problem  of  figure  3. 

The  resulting  model  is  conceptually  simple,  as  it  applies  equally  well  at  short  and  long  ranges 
(in  the  far  Geld),  unlike  either  ray  tracing  or  modal  models.  Like  ray  tracing  and  coupled 
mode  analysis,  however,  the  parabolic  equation  neglects  the  ha  tk scattered  Geld.  This  can  be 
seen  with  reference  to  equation  (12)  in  which  only  the  positive  exponential  term  eiki®  is 
included.  Other  restrictions  inherent  in  the  parabolic  approximation  are  discussed  below. 


HEIGHT  (Z) 


U(0, Z) 


RANGE  (X) 


Figure  3  Parabolic  equation  boundary  value  problem 


2J.1  Assumptions 

The  effect  of  the  paraxial  approximation  can  be  seen  by  substituting  a  trial  plane  wave 
solution  into  the  two-dimensional  wave  equation  (3,14)  and  into  the  parabolic  equation 
(13)  and  then  comparing  the  results.  If  the  wavenumber  components  in  the  x-  and  z- 
direclions  are  denoted  are  denoted  by  ^  and  k2  respectively,  then  a  plane  wave  solution 
to  equation  (14)  can  be  expressed  by 

u  -  cxpji  (k2x  +  ktz)) 

Substituting  this  into  equation  (14)  gives 

-kj1  u  +■  2ikxku  -  V’  u  +  Id  (m*  - 1)  =  0 

that  is, 

k^  *  ±  yt  (1  +  tf  -1)  -  k2* 


« 
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Now  m2  **  1  so  that 

kx-±  yt-v 

If  8  now  specifics  the  angle  of  propagation  with  respect  to  horizontal,  then  this 
equation  is  merely  a  statement  of  the  trigonometric  identity 

cos©  »  £  J  1  -  siri1  © 

The  plane  wave  solution  to  the  parabolic  equation  (equation  IS)  can  be  expressed  in 
the  form 

u  =  expfi  {(k,-  k)x  +  kj2)l 
Substituting  into  equation  (14)  gives 

-k,2  -  2(1 It,-  k)k  +  V  (of  -  1)  =  0 
which  can  be  reexpressed  as 


-k,2  -  2k,k  t-  k2  (2+m1  - 1)  =■  0 
Now  m2  ~  1  so  ,i 


which  corresponds  to 


cos© 


sin*  8 
2 


This  is  a  good  approximation  provided  that  0  is  less  than  15-20* .  This  limitation  on 
the  maximum  angle  of  propagation  is  acceptable  for  dueling  problems  since  ;l  is  only 
the  energy  at  small  angles  to  the  horizontal  (typically  smaller  than  1  degree)  that  will  be 
trapped  by  the  duct.  Therefore  energy  at  angles  greater  than  a  few  degrees  from 
horizontal  need  not  be  considered,  thereby  allowing  a  saving  in  computation  (refer  to 
section  43  for  details). 


The  second  assumption  (b)  can  be  rewritten  in  rectangular  coordinates: 


t 


_k_ 

100 


(19) 


The  term  e  /£ ,  can  be  associated  with  the  radius  of  curvature  or  propagated  energy 
rays  resulting  from  horizontal  variations  in  e .  Equation  (19)  can  then  be  represented 
in  the  form  «  /e  >  100X/2tr,  implying  (hat  the  radius  of  curvature  of  propagating 
rays  must  be  Urge  with  respect  to  a  wavelength.  In  other  words,  equation  (19)  infers 
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that  horizontal  variation  in  £  must  be  reasonably  slow.  Since  £  =  e  ^  , 
e  /e^rf  /o'  x,  which  can  be  approximated  by  rf  Sx/6  jf  .  So  (or  a  worst-case 
frequency  of  30  MHz1  and  for  a  typical  value  of  n  (say  1 .00035),  S  n?  must  be  smaller 
than,  or  equal  to,  0.00629/m,  which,  in  practice,  will  be  satisfied.  Thus  assumption  (b) 
is  valid  for  the  frequency  ranges  over  which  this  method  is  applicable.  Similarly, 
assumption  (d)  implies  that  vertical  variations  in  refractive  index  must  be  slow. 

Assumption  (c),  rewritten  in  rectangular  coordinates,  becomes 

*>5p  (2)) 

This  implies  that  the  parabolic  equation  neglects  the  near-field  effects,  since  equation 
(20)  is  only  valid  at  ranges  greater  than  sixteen  wavelengths  from  the  source. 

Like  ray  tracing  and  coupled  mode  theory  methods,  the  parabolic  equation  neglect  the 
bacicscattered  field.  It  does  however  retain  all  diffraction  effects  associated  with  tne 
propagating  medium  and  is  therefore  valid  in  regions  where  problems  with  focusing  of 
rays  and  identification  of  ray  families  cause  ray  tracing  methods  to  break  down.  The 
parabolic  equation  also  retains  the  full  coupling  between  waveguide  modes  that  mode 
theory  requires  for  horizontally  inhomogeneous  atmospheres. 


3  SPLIT  STEP  FOURIER  SOLUTION 

To  model  problems  of  tropospheric  propagation  in  ducting  environments,  a  solution  to  equation  (14) 
is  required  for  realistic  range-dependent  refractive  index  structures.  The  most  tractable  algorithm 
that  has  been  developed  to  date,  and  the  one  which  is  the  easiest  to  implement,  is  the  split  step 
Fourier  solution  developed  in  the  underwater  acoustics  arena  by  (Hardin  and  Tappcrt,  1973),  which 
uses  the  computational  speed  of  the  Fast  Fourier  Transform  to  advance  solution  over  one  range  step 
size,  fix. 

3.1  Formulation 

Equation  (15)  can  be  re-expressed  in  the  form 

«»,(**)  =  l i  (x^)  -l)  +  2^J  u(x^) 

-  i  [  A(x^)  v  B(z)  1  u(x^)  (21)  | 

1.  At  frequencies  below  VHF,  ionospheric  effects  become  increasingly  important. 
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where 


and 


A(V)  =  (w)  -  1} 


B<z)  = 


JL± 

2k  Sz> 


If  the  index  of  refraction  can  be  assumed  to  be  slowly  varying  with  range,  then  the  solution  of 
(21)  will  then  take  the  form 

u(x+4x,z)  »  i<5xexp(A+B)u  (22) 

provided  that  the  commutator 

(  (A  +  B),  |  (A  +  B)dxJ  -  0 

In  general,  this  will  only  be  true  if  A(x,z)  is  independent  of  altitude,  z  -  a  restriction  which 
would  render  this  solution  useless  for  problems  involving  ducting  structures.  Instead,  this 
condition  restricts  the  application  of  such  a  solution  to  environments  in  which  A  is  sufficiently 
slowly  varying  with  respect  to  altitude  that  the  error  incurred  by  this  assumption  is  negligible. 


The  split  step  Fourier  algorithm  involves  the  assumption  that 

u(«+4j,z)  =  t'A^IclB^,u(ipt)  (23) 

The  validity  of  this  assumption  requires  that 

[  A,  B  ]  -  0 


where  (A,B)  =  AB-BA  denotes  the  commutator  of  A  and  B.  Again  this  docs  not  hold,  since 


eA‘B  =  eA  eB  er-  eC’eC*... 


where  Ck  is  a  linear  combination  of  k-fold  commutators  of  A  and  C.  in  particular, 
c2  -  -  \  IA.B);  C3  =  i|A,[A,B])  +  |(BjA£U 


(For  discussion  of  this  exponential  identity,  refer  for  example,  to  Steinberg,  1985.) 

The  evaluation  of  exp^ 1  is  straight-forward.  Since  E'(z)  cotaains  a  second-order  differential 
operator  in  z,  the  term  clB^  "u(x,z)  can  be  evaluated  by  means  of  a  Fourier  transform 
operation  between  the  spatial  domain,  z,  and  the  spa’*?’  fieotency  variable  p  (p-k  sin  9, 
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where  0  represents  the  angle  of  propagation  from  the  horizon): 


eiB5Vv)=/'{  /fo*6'  u(x,z)l  ) 

T 

r+«o 

fMKZ)}-y~ 

n(x^z)  e'lpI  dz 

)  — CO 

"1 

r-b® 

/’K^)}  -yz 

u(jcz)  eH(B  dz 

J  —00 

So  equation  (23)  becomes 

u(*+«  *.  z)  =  o'*6  “  /'  {  / [eiBi  *  »{V)]  } 

-  e'k«»("' -I'Wft  e->p6*P  /[U(W)J  }  (24) 

Equation  (24)  is  obviously  an  iterative  solution,  where  the  field  at  any  range  can  be  evaluated 
by  advancing  the  solution  in  small  steps,  given  an  initial  field  u(0,2)  and  satisfaction  of  the 
boundary  conditions  at  the  earth's  surface  and  at  some  maximum  altitude.  The  solution  at 
each  range  is  advanced  in  two  stages:  firstly,  the  field  is  advanced  as  if  the  propagation 
medium  was  homogeneous,  thus  accounting  for  diffraction  effects  alone  (this  is  the 
contribution  of  the  forward  and  inverse  Fourier  transform  operations  of  equation  (24));  and 
secondly,  the  effects  of  the  meteorological  environment  are  applied  to  the  field  (the  tr  in  the 
first  exponential  term  of  equation  (24)  accounts  for  these).  The  vertical  refractive  index 
profiler  required  for  (his  term  are  usually  obtained  from  radiosonde  data,  or  they  can  be 
modelled  analytically.  Since  cquaiioo  (24)  is  a  marching  solution,  different  refractive  index 
profiles  can  be  entered  into  the  calculations  at  each  range  step,  thereby  facilitating  simple 
treatment  of  inhomogeneous  atmospheres.  It  is  important  to  note  that,  with  this  solution 
method,  complicated  refractive  structures  (varying  in  both  vertical  and  horizontal  directions) 
do  not  require  any  additional  computation  to  that  for  standard  atmospheres.  It  will  be 
recalled  that  for  the  other  full-wave  methods  discussed  in  section  2.2.2,  the  computational 
overhead  of  the  solutions  is  proportional  to  the  number  of  modes  required  to  model  the 
relevant  refractive  structures. 


The  error  involved  in  the  assumption  of  equation  (24)  can  be  obtained  by  differentiating  (2V 
and  comparing  with  (21).  To  third  order  terms,  the  error  is 


.  ...  dm  dmdu  mud2m 

dxl.km^u  * 


dzdz 


Zdz2 


u.dm.z. 

Zdz  ' 


<3x>*  ^(mfH1)2  ♦  0[Ox)3l  (25) 

£.  oZ 


However,  this  error  is  dependent  on  the  particular  ‘ split*  chosen  for  the  exponential  terms  in 
A  and  B.  The  regime  chosen  by  (Hardin  and  Tapper*.,  1973),  for  example,  assumes 

u(x+£  x)»  e*^ **  e*^  u(x^) 


IS 
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Toe  error  tikxjt  rr-i  by  assiiming  the  operators  A  and  B  commute  in  this  way  is 


(Sx)3,.  .am2.2  4  dV  3^u  _  4  33mZ  du  _  1  3jtiZ. 

96  *  3z  U  k  32*  32*  k  8z*  dz  k  3z‘ 


The  optimal  splitting  that  will  minimise  solution  error  is  thus  a  subject  for  further  research. 


3.2  Stability 


Id  this  section  the  requirement  of  solution  stability  is  investigated-  A  solution  perturbation 
analysis  is  included  in  section  S. 


The  solution  formulated  in  section  3.1  will  be  stable  if  the  difference  between  the  theoretical 
and  numerical  solutions  remains  bounded  as  the  solution  is  stepped  out  in  range.  From 
equation  (25)  it  is  evident  that  the  error  will  be  bounded  as  long  as  the  refractive  index 
gradients  remain  bounded.  In  other  words,  as  long  as  the  refractive  index  gradients  are  such 
as  to  accommodate  the  requirement*  of  the  assumptions  discussed  in  section  2.3.1,  the  solution 
will  be  stable. 


4  ASPECTS  OF  IMPLEMENTATION 

This  section  deals  with  implementation  of  the  solution  of  section  3.  To  date,  this  solui'on  appears  to 
be  the  most  computationally  efficient  method  that  can  be  applied  to  modelling  tropospheric 
propagation  in  ducting  conditions.  However  Other  solutions  have  been  de  , 'eloped  for  underwater 
acoustics  applications.  The  underwater  acoustics  problem  is  very  similar  to  the  tropospheric  case  but 
includes  the  complication  of  a  sea-bottom  boundary  condition.  (Kewfcy  rt  al,  1984}  found  that  the 
interaction  of  this  boundary  in  shallow  water  problems  rendered  the  split  step  Fourier  solution 
unsuitable  for  modelling  underwater  problems  of  this  type.  Several  alternative  sob -ions  have  been 
developed.  For  example,  whereas  in  I  he  sotk  step  Fourier  scheme  the  sc  coot  order  partial 
differential  operator  in  heigh:  of  the  parabolic  equation  is  represented  by  the  inverse  Fourier 
transform  of  its  Fourier  transform,  in  the  ordinary  differentia)  equation  (ODE)  approach  adopted  by 
(Lee  and  Papadakis,  1980),  this  operator  is  approximated  by  •  central  finite  difference  operator  which 
converts  the  partial  differential  equation  into  a  system  of  lint  order  ODEs  which  arc  then  solved  by 
application  of  a  nonlinear  mukistep  numerical  method.  Another  solution,  proposed  by  (Lee,  Botscas 
and  Papadakis,  1981),  uses  the  Crick-Nkoisoo  scheme  for  variable  coefficients  to  solve  the  parabolic 
equation.  These  alternate  solutions  may  prove  useful  lor  complicated  tropospheric  problems,  but  the 
increased  computational  intensity  involved  with  these,  when  compered  with  that  of  the  split  step 
Fourier  solution,  is  a  major  disincentive. 
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fhe  following  is  a  discussion  of  the  application  of  initial  and  boundary  conditions  to  the  split  step 
solution  of  section  3,  and  begins  with  a  description  of  the  way  in  which  the  surface  boundary 
condition  is  incorporated  into  the  present  implementation.  The  specification  of  the  maximum  height 
boundary  is  included  in  section  43,  "Signal  Processing  Aspects*,  as  the  particular  implementation  of 
this  condition  is  a  consequence  of  application  of  the  discrete  Fourier  transform.  At  this  maximum 
height  boundary  a  pseudo  radiation  condition  is  introduced  by  smooth  attenuation  of  the  field.  Other 
aspects  of  implementation  that  are  introduced  as  a  by-product  of  the  discrete  Fourier  transform  are 
also  included  under  section  43. 

4.1  Surface  boundary  condition 

With  the  assumption  that  the  skin  depth  of  electromagnetic  radiation  within  the  earth  is  small 
compared  to  the  earth’s  radius  of  curvature,  the  boundary  condition  at  the  surface  of  a  smooth, 
finitely  conducting  earth  can  be  approximated  by  (Senior,  1960): 

57  u(x,0)  +  ikrju(x,0)  =0  (26) 

where 

riv  -  jHt  /  (e  r  +  '&  )  (27) 

_ ^0 

0H  =7  (er+J2.)/Mr  (28) 

can  be  substituted  into  (26)  for  vertical  and  horizontal  polarisation  respectively.  In  these 
expressions  £  and  O  are  the  relative  permittivity  and  absolute  conductivity  of  the  surface 
medium  respectively,  and  £  g  is  the  free  space  permittivity,  w  is  the  radian  frequency  of  the 
propagating  electromagnetic  waves,  wbiie  t*r  refers  to  the  relative  permeability  of  the  surface 
medium,  which  is  often  assumed  to  be  unity. 

Polarisation  of  the  propagating  waves  is  incorporated  in  the  surface  boundary  condition.  From 
(27)  and  (28)  it  is  evident  that  for  a  perfectly  conducting  surface  -f-u(x,0)=0  (vertical 
polarisation)  or  u(x,0)=0  (horizontal  polarisation),  in  which  case  the  field  must  be  symmetric 
or  antisymmetric  respectively  about  the  surface.  The  surface  boundary  condition  can 
therefore  be  satisfied  using  a  linear  combination  of  even  and  odd  solutions  at  each  step,  with 
the  assumption  that  their  complex  combination  ratio  is  constant  over  a  range  interval,  Sx. 
Errors  introduced  by  this  assumption  are  negligible  provided  that  5  x  is  not  too  large. 

The  present  implementation  of  equation  (24)  splits  the  field  into  even  and  odd  parts  after  the 
inverse  Fourier  transform  has  been  evaluated,  applies  a  combination  formula  for  these  parts 
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which  satisfies  (26),  then  multiplies  both  parts  by  the  environmental  exponential  term  of  (24) 
before  reconstructing  the  full  Geld  from  the  even  and  odd  parts,  ready  to  begin  solution  at  the 
next  range  step. 

Losses  due  to  irregular  terrain  and  sea  state  conditions  can  also  be  incorporated  into  the 
model  but  these  have  not  been  addressed  in  this  report. 

42  Initial  starting  field 

Since  the  tropospheric  problem  formulated  in  section  2.1  is  an  initial  value  one,  it  is  necessary 
to  specify  the  field  at  range  x=0  in  order  to  start  the  computational  procedure.  There  are 
several  ways  that  the  starting  field  can  be  generated,  the  most  obvious  being  to  solve  the  full 
wave  equation  in  a  small  region  containing  the  source  and  to  extend  this  out  several 
wavelengths  in  range  from  the  source  to  the  region  where  the  parabolic  equation  becomes 
valid  (refer  to  restriction  e,  page  11).  If  the  troposphere  can  be  assumed  to  be  exactly 
stratified  near  the  source,  then  this  solution  can  be  obtained  by  separation  of  variables  and 
calculation  of  the  normal  modes.  This  is,  however,  a  rather  complicated  procedure 
considering  that  the  energy  of  interest  in  ducting  environments  is  merely  that  at  small  angles  of 
incidence.  Instead,  the  initial  radiated  field  can  be  calculated  with  the  aid  of  Fourier  shifting 
theory. 

From  antenna  theory,  the  radiated  field  is  proportional  to  the  Fourier  transform  of  the  field  in 
the  antenna  aperture.  Antenna  pattern  shape,  height  and  steered  direction  can  thus  be 
modelled  using  the  Fourier  shift  theorems  to  produce  the  starting  field  for  calculations.  This 
approach  has  been  implemented  in  (Dockery,  1988)  to  specify  arbitrary  symmetric  antenna 
patterns. 

The  approach  taken  by  the  author  for  example'  included  in  section  7  is  to  model  the  starting 
field  by  a  Gaussian  function  which  is  derived  to  assymptotically  match  the  elliptic  wave 
equation  Green’s  function  for  a  homogeneous  atmosphere. 

It  is  important  to  note  that  the  initial  Geld  must  be  band-limited.  The  first  reason  for  this 
arises  because  of  the  limit  on  propagation  angle  imposed  by  the  parabolic  approximation 
(refer  to  section  2J. 1  for  details).  There  is  however  a  further  constraint  imposed  by  the  fast 
Fourier  transform  algorithm,  as  discussed  in  section  43. 
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4J  Signal  processing  aspects 

In  order  to  assist  examination  of  the  signal  processing  aspects  of  the  split  step  solution,  it  is 
useful  to  rewrite  (24)  in  terms  of  dimensionless  variables: 

u(x+6  x,  z)  -  e*p{i  [£„’  (i  i)  +  i,*  (i,  i)l>  /'[expt X  )  /(“)(*  p))  (29) 

where 

z  =  fcz  p 

x  =  kx  dx  =  k$x 


are  all  measured  in  wavelengths;  and  BR2  (x,z),  H,2  (x^)  are  the  real  and  imaginary 
components  of  the  refractive  index  term  that  has  been  redefined  in  terms  of  wavelength. 

The  first  issue  of  concern  is  the  fact  that  the  solution  must  satisfy  the  boundary  condition  that 
the  field  is  attenuated  to  zero  at  infinite  height.  However,  only  a  finite  height  can  be 
accommodated  by  the  discrete  Fourier  transform.  In  order  to  avoid  unnecessarily  large  FFT 
sizes  an  absorbing  layer  is  applied  above  the  maximum  altitude  of  interest  to  slowly  attenuate 
the  field,  whilst  at  the  same  time  minimising  reflections  from  this  interface. 

The  absorbing  layer  implemented  in  the  examples  included  in  section  7  was  constructed  by 
applying  a  simple  raised  cosine  window  to  the  imaginary  part  of  the  square  of  the  refractive 
index  term.  From  the  first  ciqxinentia!  term  of  (29)  it  is  evident  that  this  will  cause  an 
exponential  attenuation.  The  raised  cosine  function  was  chosen  for  its  differentiability,  but 
other  continuous  functions  can  be  applied  for  this  purpose. 

Another  signal  processing  consideration  is  choice  of  FFT  size.  In  order  to  avoid  aliasing  in  the 
spatial  frequency  domain  (p-space),  sampling  in  height,  z,  should  be  no  coarser  than  one  half 
wavelength.  For  most  applications,  however,  this  would  require  very  large  transform  lengths. 
Since  only  energy  that  is  propagated  at  an  angle  near  to  horizontal  will  be  trapped  by  ducting 
layers,  it  is  reasonable  to  limit  the  maximum  angle  of  propagation,  8 ,  that  is  to  be  considered 
by  the  solution.  This  restriction  on©  corresponds  to  a  cut-off  in  the  spatial  frequency  domain 
(p«k  sin  8)  and  so  ignores  that  part  of  the  spectrum  which  would  Ckbcrwise  contribute  to 
aliasing.  Sampling  in  the  z-direction  can  now  be  coarser  than  the  Nyquisl  half  wavelength. 
For  a  specified  maximum  spatial  frequency,  the  maximum  height  sample  size,  6  z,  which 
would  avoid  aliasing  is  given  by  (2p(nu)  1  and  the  maximum  altitude  which  can  be  considered 
is  given  (for  a  Fourier  transform  length  of  N)  by  N£  z/2. 
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In  the  present  implementation  of  the  model,  this  low  pass  filtering  of  spatial  frequency  is 
achieved  by  applying  a  simple  raised  cosine  window  to  the  diffraction  exponential  term  of  (29) 
before  the  inverse  Fourier  transform  is  taken.  This  approach  has  been  found  by  the  author  to 
be  satisfactory  for  refractive  structures  implemented  to  date.  More  elaborate  filtering  may  be 
found  necessary,  however,  in  other  cases. 

There  Is  another  non-trivial  signal  processing  consideration.  It  is  important  to  note  that  the 
value  of  p^  used  to  avoid  aliasing  in  the  spatial  frequency  domain  should  also  take  into 
account  the  allowable  bound  in  the  spatial  domain,  z.  This  is  to  ensure  that  aliasing  does  not 
result  when  either  the  forward  or  inverse  transforms  are  taken  after  multiplication  by  the 
exponential  terms  of  (29)  in  the  spatial  and  spatial  frequency  domains.  The  required  value  of 
Pmu  therefore  be  slightly  less  than  that  derived  solely  from  N  and  5  z  as  indicated  in  the 
previous  paragraph. 


5  SOLUTION  PERTURBATION  ANALYSIS 

The  purpose  of  this  section  is  to  examine  the  consequences  of  errors  in  estimating  solution 
parameters.  The  effects  of  choice  of  range  step  size,  6  x,  and  height  sample  size,  5  z,  have  already 
been  covered  in  sections  3.2  and  4.3  respectively.  These  however  arc  parameters  that  can  be  chosen 
in  order  to  give  the  desired  accuracy  (in  the  case  of  S  x)  or  the  maximum  altitude/  angle  of 
propagation  that  is  to  be  considered  by  the  solution  (in  the  case  of  4  z).  What  is  investigated  here  is 
the  effect  that  errors  in  the  'measured*  parameters  have  on  solution.  By  'measured*  is  meant  namely 
the  refractive  index,  n,  and  initial  field,  u(0,z),  components;  these  arc  the  inputs  to  the  system. 

£.1  Perturbations  in  refractive  index 

Let  the  modified  refractive  index  term  of  (24)  be  perturbed  by  a  complex  amount  a  and  let  the 
corresponding  solution  at  range  x+S  x  be  denoted  by  u’n  j.  Equation  (24)  can  then  be  written 
for  the  perturbed  and  unperturbed  cases  as 

Vi  =  'ikra'4l/J  /*< cV  Sx/1  A“„>  > 
u’nM  =  c«dn>'-a*«/2  /»{  ev  $x/2  /(u-n)  , 

where  m’-m’  -I.  Then 

=  *««'***'*  /‘I  *V  6 'n I /<«•.>/(“.>  I  > 

+  (e*f«’«<*>*x/2.e*m*i/i)  eV  <5 1/2  j 
“  v,  +  v2  (29) 
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where  v,,  v2  correspond  to  the  terms  of  this  sum.  Now 

Wjl  =  e'k{te(m'*a,,i,/2  |/‘{  eV  Sx/2[f(*'j-faj])  | 

where  (m(  )  denotes  the  imaginary  part  of  (  ).  So,  finding  the  norm2  of  v}  with  respect  to 

the  measure  jjp , 

yar 

l|v,||  5  nu’.-u.H 


Similarly,  the  Lj  norm  of  is 

IIVjll  S  ■?{  le^'6*/2!  |e“4^-l|  }  llu.ll 
<  SL*  “]a|  ||Uii|| 

Now,  by  applying  the  triangular  inequality  to  (29), 

HU’»+l-Un*lll  -  "Vt'i  +  l  'V2«  I 

<  e-k[mm(Im(ni'+a)))6l/2  |  |u>  .  u  |j 

+  e-M»«(l»(«')})dr/J  ««|a|  ||n^(  (30) 

Let  e  n=|  |u’a-u0l|,  e  1“  llu,B»1-uB+lll-  Then  equation  (30)  is  simply  a  difference  equation 
of  the  form 


Recursively  expanding  this  equation  shows  that 

Now 


A  <  e"<™"ll<»(*')})4  «/*  e-kJmM{Im(a)}|4«/2  <  j 


and 

Bb  =  »/*  &  ™  |Qj  ||Ua|| 


2.  For  discussion  of  the  continuity  properties  of  the  Fourier  Transform  consult  (Desoer  and  t 

Vidyasagar,  1S75). 
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Since  ||ul|  is  bounded,  B  will  be  bounded  for  all  So  fot  A,  <1,  €  ->0,  and  for  A.  =  1,  it 

OB  Bun 

is  required  that  £  “  |a  (aS  t^l)  |  <®. 

A 

This  implies  Chat  small  perturbations  in  refractive  index  will  tend  to  zero  as  the  solution  is 
stepped  out  in  range  and  will  not  cause  the  solution  to  become  unstable. 

52  Perturbations  in  initial  field 

It  is  also  important  to  consider  perturbations  in  the  initial  field  since  these  correspond  to 
errors  in  calculation  of  the  antenna  field  or  to  inaccuracies  in  the  antenna  pattern  entered  into 
calculations. 

Rewriting  equation  (24) 

Vi  = 

V,l  =  |/‘{  e1^  il/2(/(uii)l  }  | 

The  Lj  norm  is  then 


I IV ,11  5  l|uDH 


Recursively  expanding  this  equation  shows  that 

MU, 1 1  <  {  V"?  j-MlnKm1  ||U()|| 


and  so 


lluoll  S  e4*^2  HUgll 

where 


q  =  (w)) 

So  the  solution  will  be  stable  if  0  and  assymptotically  stable  if  q<0.  q  will  only  be  positive  if 
energy  is  added  to  the  system  which  is  not  the  case  for  the  dissipative  nature  of 
electromagnetic  waves  in  the  troposphere.  Hence  perturbations  in  the  initial  field  will  not 
result  in  instability. 
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6  REFRACTIVE  VARIATIONS  CAUSING  DUCTING 

Between  3  kHz  and  3000  ghz  'be  index  of  refraction  is  essentially  independent  of  frequency.  This 
index,  n,  is  a  value  very  dose  to  unity  and  another  index  -  radio  refractivity,  N  -  is  commonly  used  for 
mathematical  convenience: 

N  -  (n  -  1)  106 

-  77.6  P/T  +  3.73x10s  e/T  (COR,  1986) 


where  P  -  atmospheric  pressure  (mb) 

T  =  absolute  temperature  (K) 
e  =  water  vapour  pressure  (mb) 

In  'standard'  atmospheric  conditions  each  of  these  parameters  (atmospheric  pressure,  temperature 
and  water  vapour  pressure)  decreases  with  increasing  altitude.  However,  meteorological  conditions 
often  deviate  from  this  standard  atmosphere  and  an  inversion  layer  of  increasing  refractivity  with 
altitude  may  form.  Radio  refractivity  gradient  varies  primarily  due  to  variations  in  temperature  and 
water  vapour  concentration,  such  variations  occurring  oo  time  scales  ranging  from  seconds  to  months. 
Values  of  these  quantities  are  usually  derived  from  radiosonde  measurements,  from  which  vertical 
profiles  of  refractive  modulus,  M,  can  be  constructed.  M  is  defined  as 

M  =  N  +  IC^b/a 

where  h  is  height  above  the  earth’s  surface  and  a  denotes  the  earth’s  radius.  Refractive  modulus 
includes  both  atmospheric  refraction  and  a  correction  term  to  account  for  the  effects  of  the  earth’s 
curvature. 

When  the  gradient  of  refractive  modulus  with  respect  to  height  is  negative,  ducting  is  present. 
Figures  4-6  contain  three  idealised  refractive  modulus  profiles,  each  representing  different  categories 
of  ducting.  Figure  4(a)  illustrates  the  presence  of  an  evaporative  duct,  formed  by  an  inversion  layer  at 
the  surface  of  the  earth.  The  duct  extends  from  the  earth  surface  to  the  top  of  the  inversion  layer,  as 
indicated.  This  type  of  duct  is  regularly  found  over  warm  bodies  of  water  and  ts  generally  caused  by  a 
temperature  inversion  near  the  surface  that  is  accentuated  by  intense  relative  humidity  as  a  result  of 
evaporation.  When  this  type  of  duct  occurs  over  land  ■>  is  usually  referred  to  as  a  surface  duct.  As 
figure  4(b)  illustrates,  energy  propagated  from  within  the  duct  at  a  few  degrees  from  horizontal  will 
be  trapped.  Evaporation  duds  axe  typically  less  than  10  m  thick  and  surface  values  of  refractive 
modulus  around  the  world  range  between  300-400  M  units. 
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Figure  5(a)  illustrates  a  profile  in  which  an  elevated  duct  is  present.  Such  profiles  contain  two 
inflection  points  above  the  surface,  each  with  a  refractive  modulus  value  greater  than  that  at  the 
surface.  The  duct,  indicated  on  the  profile  by  a  vertical  line,  extends  from  the  top  of  the  inversion 
layer  to  the  height  at  which  the  refractive  modulus  equals  that  at  the  upper  duct  boundary.  As 
illustrated  by  figure  5(b),  only  energy  launched  into  the  duct  at  shallow  angles  will  be  trapped. 
Elevated  ducts  range  in  thickness  bom  a  few  metres  (these  tend  to  affect  propagation  above 
microwave  frequencies)  up  to  hundreds  of  metres  (affecting  propagation  at  frequencies  above  YHF). 

Figure  6(a)  shows  a  surface-based  elevated  duct  The  duct  is  elevated  since  two  inflection  points  are 
present  in  the  profile,  but  it  is  surface-based  since  the  refractive  modulus  at  the  upper  inflection  point 
is  less  than  that  at  the  surface.  Energy  is  trapped  by  the  duct  as  indicated  by  the  ray  paths  of  figure 
6(b). 


The  ray  traces  in  figures  4-6  give  an  indication  of  the  behaviour  of  energy  in  and  around  ducts,  but  do 
not  indicate  the  relative  strengths  of  that  energy.  Ducts  which  include  both  transmitter  and  receiver 
(target)  will  generally  enhance  the  received  signal,  while  the  duct  may  act  to  shield  the  signal  if 
transmitter  or  receiver  (target)  is  situated  outside  the  duct.  Since  refractive  index  commonly  varies 
with  range  as  well  as  with  altitude,  quite  complex  ducting  structures  are  not  uncommon,  from  which 
complicated  field  patterns  emerge.  What  would  be,  at  one  specific  altitude  and  range,  an  optimal 
location  from  which  to  communicate  with  a  specified  receiving  point,  could  at  a  nearby  position  be 
quite  shielded  from  the  receiver.  The  advantage  to  be  gained  by  a  communications  or  surveillance 
system  operator  through  a  reliable  propagation  prediction  technique  is  therefore  obvious.  The 
following  examples  of  calculated  field  patterns  for  various  ducting  structures  serve  to  illustrate  the 
complicated  nature  of  '.he  effects  on  propagation  that  ducting  can  cause. 


Figure  4  Ray  paths  and  stylised  vertical  profile  for  evaporative  duct 
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FigurcS  Rujr  paths  nod  stylised  vertical  profile  for  derated  duct 


u 

u>  <w 

Figure  6  Kay  paths  and  stylised  vertical  profile  for  surface-based  elevated  duct 


7  COVERAGE  DIAGRAM  EXAMPLES 

The  split  step  Fourier  solution  to  the  parabolic  wave  equation  has  been  implemented  by  the  author 
on  a  personal  computer  with  an  in-board  fast  parallel  processing  card.  The  program  is  written  in 
FORTRAN  and  uses  a  graphics  display  package  called  HGRAPH  to  produce  coloured  coverage 
diagrams  of  relative  field  strength  on  a  height  versus  range  flattened -earth  grid.  Execution  times  are 
directly  proportional  to  the  range  out  to  which  calculations  are  extended,  and  depend  on  the  FFT 
sizes  in  the  solution  which  are  in  turn  dependent  oo  frequency.  As  an  indication,  the  examples  In  this 
section  for  the  1  ghz  source  took  approximately  two  hours  to  generate  for  a  range  of  400  km  (FFT 
size  was  3192  points),  while  lor  the  1UU  MHz  source,  plots  took  approximately  twenty  minutes  over 
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the  same  range  (with  FFT  sizes  of  1024  points)-  As  indicated  in  scctioo  43,  the  discretisation  size  in 
altitude  (sod  therefore  the  length  of  the  FFT)  is  related  to  the  maxim  am  allowable  spatial  frequency 
in  p- space.  This  maximiun  spatial  frequency  directly  relates  to  die  maximum  angle  of  propagation 
considered  by  the  solution.  So  if  energy  at  angles  that  are  greater  than  the  marimnm  is  to  be 
considered,  then  the  altitude  discretisation  size  win  need  to  be  decreased  and  therefore  larger  FFT 
sizes  will  be  required.  This  will  obviously  increase  execution  times. 

Note  that,  unlike  for  other  full-wave  methods,  execution  time  is  effectively  independent  of  duct 
profile  complexity.  The  qualifying  'effectively'  here  is  included  because  there  may  be  cases  in  which 
the  duct  profiles  need  various  amounts  of  interpolation  to  smooth  discontinuities,  in  the  duct 
structures  so  as  to  reduce  reflections  from  such  irregularities.  Such  increases  in  execution  time, 
however,  win  be  independent  of  frequency  and  FFT  lengths. 

The  examples  included  in  section  7.1  all  deal  with  range-independent  ducting  problems.  That  is,  the 
ducting  profiles  used  in  these  examples  are  horizontally  homogeneous  (they  vary  in  height  but  not  in 
range).  These  diagrams  illustrate  the  effects  on  forward  coverage  that  surface  and  elevated  ducts 
cause.  In  section  7 2,  these  examples  are  extended  to  include  problems  that  are  range-dependent. 
Simple  analytical  structures  are  used  to  illustrate  the  complicated  nature  of  the  patterns  resulting 
from  such  varying  ducts. 

7.1  Range-independent  problems 

Figure  7  represents  the  coverage  of  a  1  gfrz  vertically  polarised  Gaussian  beam  source  at  a 
height  of  ISO  m  in  a  "standard"  exponential  atmosphere  (as  defined  by  CC1R  report  563-3, 
1986).  A  perfectly  conducting  earth  surface  is  assumed  and  the  subsequent  interference  lobing 
is  apparent.  Figure  8  depicts  this  same  source  in  a  horizontally  homogeneous  elevated  duct. 
The  duct  extends  from  approximately  50  m  to  600  m  and  is  of  strength  55  M-units,  as 
illustrated  in  the  vertical  refractive  index  profile  beside  the  coverage  diagram.  The  energy 
trapped  by  this  duct  is  apparent.  With  reference  to  figure  8,  it  is  dear  that  coverage  is 
enhanced  for  heights  below  600  m.  Above  this  altitude,  however,  coverage  is  significantly 
reduced  within  200  km  from  the  source.  At  larger  ranges,  there  is  some  leakage  from  the 
upper  boundary  of  the  duct.  Note  that  unlike  ray  tracing  techniques  which  tend  to  pass  or 
reflect  all  of  the  energy  at  a  duct  interface,  this  parabolic  approach  permits  partial  energy 
reflection. 

Another  point  which  is  important  to  note  with  reference  to  figure  8  is  that,  as  a  receiver  (or 
target)  moves  through  the  duct,  the  received  field  varies  as  a  Function  of  height  as  well  as  of 
range.  This  is  because  of  interference  between  waves  which  have  undergone  different 
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amounts  of  refraction/reflectioo,  and  also  (to  adopt  ray  tracing  terminology)  because  of 
focusing  effects.  This  point  has  been  graphically  illustrated  in  (Slingsby,  1990)  where,  for  a 
horizontally  homogeneous  environment,  radar  signal  returns  are  calculated  using  the 
modelling  technique  described  in  this  report.  The  resulting  plots  of  relative  returned  signal 
versus  antenna  height,  and  of  relative  returned  signal  versus  target  height,  highlight  the 
complicated  effects  that  ducts  have  on  signals  and  emphasise  the  need  for  an  accurate 
prediction  model.  Thus,  the  knowledge  that  a  duct  is  present  is  not  sufficient  to  enable  an 
operator  to  determine  at  which  position  a  transmitter  or  receiver  should  be  placed  in  order  to 
optimise  system  performance. 

In  figure  9  the  homogeneous  ducting  layer  of  figure  8  has  been  reduced  in  strength.  The 
vertical  refractive  index  profile  for  this  diagram  shows  that  this  action  has  the  effect  of  raising 
the  hase  height  of  the  duct  and  this  is  reflected  in  the  coverage  diagram.  Because  the  duct  is 
weaker,  more  energy  escapes  from  the  duct,  as  illustrated. 
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Figure  10  Coverage  diagram  for  I  GHi  source  la  inhomogeatous  duct 
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7  2  Range-dependent  problem* 

Figure  10  depict*  the  coverage  for  this  same  1  GHz  source,  only  this  time  in  a  non- 
horn  ogeneo  us  range-dependent  environment.  Study  of  this  coverage  diagram  in  light  of  the 
coverage  diagram  of  figure  9,  provides  an  indication  of  the  inaccuracy  incurred  by  assuming 
horizontal  homogeneity.  Both  diagrams  begin,  at  range  x=0,  with  an  elevated  duct  of  20  M- 
units  at  600  m.  Whilst  the  dnet  of  figure  9  remains  constant  for  all  ranges,  the  duct  of  figure  10 
has  been  weakened  at  a  rate  of  3  M -units  per  100  km.  Now  this  refractive  modulus  gradient  is 
very  slow;  experimental  results  obtained  for  elevated  ducts  over  the  Timor  Sea  by  (Barton, 
1973),  for  example,  indicate  gradients  of  32  M-nnits  per  100  km.  However  the  differences  in 
coverage  patterns  are  evident.  With  increasing  range,  much  more  energy  begins  to  escape 
from  the  top  of  the  duct  of  figure  10.  As  a  result,  the  signal  detected  by  a  receiver  placed 
350km  from  the  source  at  an  altitude  of  800  m,  is  20  db  stronger  for  the  inhomogeneous 
ducting  environment  of  figure  10  than  for  the  case  depicted  in  figure  9.  But,  on  the  other 
hand,  a  receiver  400  km  from  the  source  at  an  altitude  of  630  km  would  detect  a  signal 
approximately  20  db  weaker  than  would  be  predicted  using  the  homogeneous  ducting 
environment  of  figure  9. 

The  example  given  above  highlights  the  fact  that  methods  which  use  ’typical’3  or  ’interpolated’4 
vertical  refractive  index  profiles  to  represent  inhomogeneous  ducting  environments  are 
susceptible  to  large  errors  in  predicted  signal  coverage.  The  ability  to  model  two- 
dimcnsionally  inhomogeneous  structures  is  clearly  important.  Making  the  assumption  (hat  a 
duct  is  horizontally  uniform  can  cause  significant  underestimation  of  the  level  of  leakage  from 
the  duct  and  can  lead  to  field  strength  predictions  that  are  substantially  in  error.  Therefore, 
when  elevated  ducting  structures  are  present,  the  most  accurate  model  is  the  one  prescribed  in 
this  report. 


8  CONCLUSION 


This  report  demonstrated  the  application  of  the  parabolic  equation  method  to  range-independent  and 
range-dependent  tropospheric  propagation  modelling  problems.  The  parabolic  equation  is  an 
approximation  to  the  Helmholtz  wave  equation  which  allows  progressive  calculation  of  the 


3.  This  refers  either  to  a  vertical  refractive  index  profile  derived  from  measurements  made  at  a  single 
point  along  the  propagation  path  or  to  a  profile  that  is  estimated  by  averaging  profiles  at  several 
points  along  the  path. 

4.  Commonly,  a  linear  interpolation  between  vertical  refractive  index  profiles  is  used  to  determine  a 
representative  profile  for  a  particular  range.  Various  authors  perform  what  can  appear  to  be 
something  of  a  ’fudge’  as  they  apply  different  methods  of  interpolation  to  these  profiles  in  an  attempt 
to  verify  experimental  propagation  results. 
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propagated  electromagnetic  field  as  soiutioo  is  stepped  out  is  range.  The  particular  solution  method 
described  in  this  report  and  implemented  at  ERL  derives  from  a  technique  developed  by  (Hardin  and 
Tappert,  1973)  for  application  to  underwater  acoustics  problems.  This  ’split  step  Fourier’  soiutioo 
involves  the  discretisation  of  the  field  with  respect  to  height  and  the  use  of  the  Fast  Fourier  transform 
at  successive  range  points.  The  solution  is  stable,  with  errors  being  dependent  on  wavenumber,  range 
step  size  and  on  the  gradients  of  refractive  index  with  respect  to  height  and  range.  It  has  been  shown 
in  this  report  that  the  solution  also  behaves  well  for  perturbations  applied  to  the  initial  field  and  for 
small-scale  discontinuities  in  refractive  index  profiles. 


The  propagation  modelling  method  presented  in  this  report  has  a  number  of  advantages  over  other 
methods  used  to  evaluate  the  effects  of  ducting  layers  on  tropospheric  propagation.  The  parabolic 
equation  method  retains  all  the  diffraction  effects  associated  with  the  propagating  medium  and 
therefore  is  valid  in  those  regions  where  ray  tracing  techniques  break  down.  Unlike  ray  tracing 
methods  which  give  a  qualitative  picture  of  the  effects  of  ducting  on  coverage,  the  method  presented 
in  this  report  uses  a  full-wave  approach  to  calculate  the  amplitude  function  of  the  propagating  signal. 
The  method  differs  from  coupled  mode  techniques  in  that  it  readily  accommodates  atmospheres 
which  are  horizontally  inhomogeneous.  Whereas  for  coupled  mode  techniques  computational 
overhead  is  proportional  to  the  intricacy  of  the  duct  profiles,  computational  intensity  of  the  technique 
described  in  this  report  is  independent  of  the  duct  profile  complexity.  By  inviting  a  'marching' 
solution  the  parabolic  equation  offers  large  computational  saving  over  existing  coupled  mode 
techniques. 

The  examples  given  in  section  7  of  coverage  diagrams  generated  with  the  parabolic  equation 
technique  serve  to  illustrate  the  complicated  nature  of  the  propagation  patterns  which  arise  due  to 
ducting  layers.  It  is  particularly  important  to  note  the  difference  in  coverage  generated  for  a  source 
in  a  horizontally  inhomogeneous  ducting  atmosphere  when  compared  with  that  same  source  in  a 
homogeneous  ducting  atmosphere.  Examples  of  this  type  highlight  the  realistic  need  for  a 
mathematical  model  such  as  the  one  described  in  this  report  to  accurately  illustrate  the  effects  that 
ducting  layers  have  on  propagation  performance. 

Whilst  simple  comparisons  between  the  parabolic  equation  method  and  other  existing  techniques 
have  been  made  by  the  author,  verification  against  experimentally  obtained  propagation  data  has  not 
yet  been  achieved.  The  reason  for  this  is  that  very  little  experimental  data  is  in  existence.  The 
atmospheric  data  necessary  for  the  refractive  index  profiles  that  arc  input  into  the  model  is  usually 
collected  via  radiosonde  releases  from  the  ground  or  via  refractometer  measurements  from  an 
aircraft  .  This  can  be  a  very  expensive  process  and  to  date  there  exists  very  little  data  of  this  type  that 
is  concurrent  with  propagation  measurements.  This  point  raises  the  question  of  bow  much  practical 
use  such  a  mathematical  model  will  have  given  present  difficulties  in  attaining  the  input  data.  The 
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answer  lies  in  the  feasibility  of  remotely  sensing  this  data  by  some  other  means  than  is  currently 
employed  One  possibility  is  the  derivation  of  these  profiles  from  meteorological  satellite  data.  If 
this  remote  sensing  can  be  achieved  then  the  modelling  technique  described  in  this  report  will  allow 
enhanced  real-time  evaluation  of  communications,  surveillance  or  electronic  warfare  system 
performance. 
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